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We use a set of Schwinger-Dyson equations for the Ising Model to check several random number 
generators. For the model in two and three dimensions, it is shown that the equations are sensitive 
. . . tests of bias originated by the random numbers. The method is almost costless in computer time 

00 , when added to any simulation. 
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I. INTRODUCTION 

Among the sources of systematic error in Monte Carlo (MC) simulations, the most frightening is the lack of 
randomness in the used pseudo random number generator (PRNG). Indeed, in a modern MC simulation as much as 
^ ■ 10"'^^ random numbers may be generated Q. This is as long as the longest test of random numbers we have heard 
0^ " of 1^. Therefore, a PRNG needs to be fast and thus not too sophisticated, but it also should not bias the simulation 
results. Shift-register PRNG's |^| have become very popular, due to their speed, but they have been shown to be 
unreliable for some applications |J. The study of the trustworthiness of a PRNG is quite difficult as the answer is 
problem-dependent, algorithm-dependent and (most important) precision-dependent. For instance, in ref. some 
OO , commonly used shift-register PRNG's were shown to yield incorrect results for the two dimensional Ising model 
■ simulated with the Wolff's single-cluster algorithm Of course, this failure is related to one's statistical accuracy 
+^ ' (all the generators in ref. would be "correct" with 5% errors). In particular, the R250 shift-register PRNG was 

a found to be very dangerous for single-cluster update, but safe for its use with the Metropolis algorithm. Not long after 
[ that, R250 was shown to fail in the Metropolis update of the Blume-Capel model for some lattice sizes |^]. Another 
• example of the difficulty in certifying PRNG's can be found in ref. Q. There, RANF (the standard Cray PRNG) 
^ is shown to be "very good" in the author own wording. This means that the longest carried run did not find bias in 
O ' a two dimensional Ising model simulation, where comparison with the exact solution is possible j^. Nevertheless, it 
has produced awfully wrong results in a U(l) lattice gauge-theory simulation ||l0|] . Moreover, it is fairly common that 
one's simulation is, itself, the longest run ever carried for this particular problem (otherwise, why bother doing it?). 
Unless independent, algorithmically different simulations were performed, it is clear that one's result will not be yet 
established. Further confidence can be obtained if sensible consistency tests are also carried. In this paper we want 
I to show that Schwinger-Dyson identities may be useful in this respect, specially when no exact solution is at hand. 
Let us finally mention that the investigation of the reasons for PRNG induced bias is interesting in itself , but it 
has not yet reached predictive power (one wants to know before carrying the simulation). 



II. THE EQUATIONS 

Generally speaking, Schwinger-Dyson equations are relations of the type 

where O is an arbitrary operator and H is the hamiltonian (notice however that for Eq. (^) not to be a trivial 
= statement, O should be an odd operator if H is symmetric under the (j) —> —cf) transformation). The problem 
is that the longest MC runs are usually done in discrete spin models, for which there are no continuous variables. 
Nevertheless, for spin models the measure usually has a Z2 symmetry, which allows to obtain equations analogous 
to (^). As an example, let us consider the Ising model on the cubic lattice, with nearest neighbors interaction. The 
Hamiltonian is 
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where a are the usual Z2 spin variables. Let us call 5^ to the sum of the spins coupled to spin (7^. The self-evident 
relation 

cr=-14 (T=-l,l 

yields for any observable depending on the spin ai (and possibly also on others), 0{ai \ ...), the following relation 

(O(a,;...)) = (O(-a.;---)e-''''^'^0- (3) 

In particular, one gets 

l = (e-2/3--S,^ ^ (4) 
{a,aj) =-(a,a,e-2/3-.s.) + 25,,, (5) 

where is the Kronecker symbol. In order to gain statistics, it is useful to sum Eq. for all the lattice sites (the 
lattice size being L, its volume \sV — L^)- One obtains: 

l.^lEe-^-^^y (6) 

Summing to the nearest neighbors in Eq. (|^), we obtain an expression which has been very useful in MC Rcnor- 
malization Group investigations of the dynamics of the Poliakov loop in lattice gauge-theories : 

o = (^E-^^^(i+^"''"''))- (7) 

It is trivial to generalize Eq. (0) when more couplings are included in the Hamiltonian, as needed in a MC Renormal- 
ization Group study. 

In addition, a non-local identity is obtained from Eq. (0) summing to all i and j: 
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At this point, it is natural to ask if the right-hand side of Eqs. (j^jQ® '^^'^ be measured with reasonable statistical 
accuracy. We shall see that the answer is positive. Then the next natural question to ask is if a PRNG inducing bias 
also spoils the fulfillment of these equations. We shall find a positive answer only for Eqs. (|^) and (^. 

Finally, let us mention that the Z2 symmetry is embedded in the symmetry of many other models, therefore 
Eqs.(p|j7P) hold as they are for O(A^) spin-models, or, with trivial modifications, for SU(2iV) lattice gauge-theories. 



III. NUMERICAL RESULTS 



We have studied the Ising model (with periodic boundary conditions) in two and three dimensions at their critical 
points. Three update methods have been considered: Metropolis the Swendsen-Wang cluster method and 
Wolff's single-cluster (SC) [||. For each update, we have employed three PRNG. One has been the problematic ||] 
R250 

^R250 ^ (^R250^ + ^R.50^) ^^^j (9) 

The second has been the Parisi-Rapuano (PR) PRNG Q , which has been found not quite correct in four dimensional 
site-percolation ||l|: 

X™ = r„ xory„_6i, (10) 
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where 



= {Yn-24 + i"n-55) mod 2^^. 

Our last generator is defined with the help of a congruential generator: 

Zn+i = (16807 Z„) mod (2^^ - 1). 

Then, the PRC PRNG § is defined as 

= + 2Z„)mod23^ (11) 

Our statistics have been the following. In two dimensions we have considered a L = 16 lattice. We have simulated at 
the exact critical point up to 6 digits 

Pc = 0.440687. 

We have measured every 20 Metropolis sweeps or 20 single-clusters, performing 8 x 10^ Metropolis full-lattice sweeps, 
and updating 4 x 10^ clusters. For the Swendsen-Wang algorithm, we measure every 5 sweeps, and have generated 
the clusters 4 x 10^ times. 

In three dimensions, the critical coupling is known with great accuracy [^. We have simulated at 

P = 0.221654. 

As shown in ref. |Q, it might happen that the bias only appears for some lattice sizes. Therefore, we have studied 
i = 16 and 24 lattices. For Metropolis or single-cluster, we measure every 10 sweeps. We perform 10^ Metropolis 
sweeps, and generate 10^ clusters. In the Swendsen-Wang case, we measure every 4 sweeps, generating the clusters 
4 X 10^ times. We have found quite clear results for the different simulations, except for the single-cluster update of 
the 16^ and 16^ lattices, with PR as PRNG. We have found convenient to extend these two simulations, although 
this is in principle a dangerous procedure. Of course, one cannot proceed the run until the results "looke nice" , since 
this would bias the results. To avoid subjective decisions, we have fixed a priori the total (much longer than the 
initial) simulation time: these two simulations have been 40 times longer than the others. In this way, error bars 
shrink enough to distinguish between a large statistical fluctuation and a systematic error. 

Before presenting our results, a word of caution is in order. We have carried 27 independent simulations (3 lattices 
X 3 PRNG X 3 updates), so, the number of expected data points which are more than one standard deviation away is 
uncomfortably large. Specifically, one can easily estimate that the number of points that are more than 1.7 deviations 
away (10% probability) must be between 2 and 4. Moreover, errors are not obtained with perfect accuracy. Allowing 
a 10% error in the error determination, we have considered deviations larger than 3.3 error bars as a significant signal 
of bias (less than a 0.13% probability). 

Let us first discuss our results in two dimensions. In the left-hand side of figure |l| we plot the deviations of the 
energy and the specific- heat from their exact values We find significant deviations only for the single-cluster 
update when using R250 and PR as PRNG's (the former is not surprising |^). It is clear that the exact solution is 
the best of possible tests, but we would like to confront it with the Schwinger-Dyson test. For this, let us define the 
quantities: 




which are the right-hand side of Eqs.(^,0). Unfortunately Eq. (^) has been found to hold within errors in all cases. 
In the above expressions ( )mc is the MC average, not the expectation value. We show our results for Ai and A2 in 
the right-side of figure |l|. The only significant deviation found is in the single-cluster update with R250. This docs 
not mean that the Schwinger-Dyson identities can be fulfilled with a biasing PRNG, as this is of course a matter of 
accuracy. In fact, performing a 40 times longer run with PR, we find 

Ai = 1.00024(4), 
A2 = -0.00047(8). 
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Thus, both the exact solution test and the Schwinger-Dyson identities test are failed by this SC-PR combination, 
but the exact solution test is more sensitive in this case. 

We can discuss our results more quantitatively. For small bias, it is natural to expect that its main effect can 
be described as a shift on the coupling, from (3 to (3' = (3 — A/3. With this assumption, we can relate the different 
bias. Let AO be the the difference between the mean value of O obtained with some MC simulation, and its true 
Boltzmann average, we obtain to first order in A/3 |^ 

AO«-^AA,. (13) 

In this way, we can understand that the bias for the energy has opposite sign that the one for Ai (see footnote), 
and it is also opposite to the bias for the specific heat (it is well known that the maximum of the specific heat of 
the two dimensional Ising model in a finite lattice is at /? < /3c). The only evidence that we can offer for Eq. (|l|) 
is empirical, and it is shown in table |. Nevertheless, we find the agreement quite satisfactory for such a rough 
calculation. Moreover, from table | we can estimate that 

Ai? « -1.6 AAi, 

AC„ « 40 AAi, ^ ' 



where the coefficient for the energy is really 1.58(19), to be compared with 1.33 from Eq. (|13|). Notice that if Eq. (|1J 
could be rigourously established, it would be enough to estimate the failure in the Schwinger-Dyson test for any 
PRNG, to get the safe accuracy level for every observable. However, to our knowledge, such an interesting property 
has not been proved for any PRNG test. 

For the three dimensional case, we plot our results for the energy and the specific-heat in the left-hand side of 
figures H and |[ In this case, we unfortunately lack an exact solution to control for the bias. However, we can study 
the statistical compatibility of our data. From the plot it is apparent that the SC-R250 results are biased. In the 
figures we show the data with a weighted estimate of the energy and the specific- heat (excluding the SC-R250 data). 
In fact, no further significant deviations are found. For the Schwinger-Dyson test, we find again strong signal of bias 
in Ai and A2 for the combination of R250 with single-cluster. We also find worrying deviations in the single-cluster 
update with PR as PRNG for L — 16. To clarify if a bias is present in this case, we have performed a 40 times longer 
run. The new results are 

Ai = 0.999966(8), 
A2 = 0.000101(26). 

Thus, the PR PRNG does produce biased results in combination with the single-cluster update. In this case, we 
cannot check Eq. (|l^) directly, as the deviations in the energy and specific-heat for SC-PR are not large compared to 
the errors. Nevertheless, we can compare the bias for the SC-R250 in the 16'^ and 24^ lattices (see table ||), which is 
a test of the L dependence of the linear coefficient in Eq. (|l^). From Finite-Size Scaling theory, we can estimate that 



(Aa/Ai;)L=24 /24 



1/^ 

wl.9... (15) 



where v « 0.63 is the critical exponent for the correlation-length. From table ||, the above quotient can be estimated to 
be 1.7(5), which is certainly compatible with our prediction, but the error is so big that this is not a compelling evidence 
for Eq. (^3|). Now, if we assume again Eq. (p^), we obtain for the 16^ lattice l\E « -5AAi and AC^ w 194AAi 
From these relations, and from the estimate of AA^'^~^^ we obtain for the bias (the statistical errors in fig. |^ being 
ue and ac) 



A£;sc-PR ^ 0.00015 , cTfi = 0.00019 



AC^c-PR ^ 0.007 , ac = 0.005. ^^^^ 



Thus, it is not surprising that the bias does not show up in Fig. |[ For the 24^^ lattice we lack an accurate measure 
of the bias for Ai, and so we cannot obtain a bias estimate. 



^To avoid confusions let us recall that in this expression = ( 2^7 X^i '^i'^i}, and it is a growing function of /?. 
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As a final remark, notice that tiie sign of the bias seems to be independent of the lattice size and the space dimension 
for R250. This seems to be consistent with the simple (unidimensional) model proposed in ref. [0. However in the 
PR case, the (much smaller) bias changes sign when going from 2 to 3 dimensions. This suggests that the reason for 
bias is more involved in this case. 

IV. CONCLUSIONS 

In this work, we have shown that some Schwinger-Dyson identities, Eqs. (HJ^), are sensitive test of PRNG induced 
bias. Most important, they can be used when no exact solution is at hand. We have provided some empirical 
evidence for a simple relation between the bias induced in the different observables (our Eq. (p^)). This relation is 
obtained under the assumption that the main effect of the bias is to produce a shift on the coupling. It might be 
possible to justify this in terms of relevant and irrelevant operators, in the framework of the Renormalization Group. 
Furthermore, this suggests that an investigation along the lines of Ref. could be useful to establish which new 
couplings are generated by the PRNG-induced bias. If this relation could be established, the SD equation test would 
provide an estimate on the maximum safe accuracy that one can get for any observable, with the given PRNG. 

In three dimensions, where there is no exact solution at hand, the Schwinger-Dyson Equations test has shown that 
the single-cluster update with the R250 and PR PRNG's produces biased results, without resource to seven more 
simulations. It should be noticed that the measure of the Schwinger-Dyson equations comes almost for free, as the 
number of possible exponential factors is finite, and the local energy should be measured anyway. Disk storage is not 
a shortcoming either, because no reweighting is to be done, and the calculation can be made "on the fly" . They 
are also extremely helpful for code debugging. So, we believe Schwinger-Dyson equations to be very useful tools, 
which can be easily measured in almost every circumstances. 
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TABLE I. The bias for £, C^,, Ai, A2, for the 16^ lattice simulated with the SC-R250, SC-PR combinations. To be able 
of measuring the bias in the SC-PR combination we have needed a much longer simulation (see text). The constancy of the 
ratios is a check for Eq. (|l|). 





AE 


AC, 


AAi 


AA2 


SC-R250 


-0.00235(11) 


0.0599(14) 


0.00148(19) 


-0.0029(4) 


SC-PR 


-0.00057(2) 


0.0115(2) 


0.00024(4) 


-0.00047(8) 



Ratio 0.242(14) 0.192(5) 0.16(3) 0.16(4) 



TABLE IL Bias for E, Ci,, Ai, A2, for the 16^ and 24^ lattices simulated with the SC-R250 combination. The "correct" value 
has been taken from the averaged estimate of Figs. ^ and ^ The ratios test the lattice-size dependence of the coefficients in 
Eg. 



L 


AE 


AC, 


AAi 


AA2 


16^ 


-0.00124(19) 


0.051(4) 


0.00026(5) 


-0.00082(16) 


24^ 


-0.00066(13) 


0.046(5) 


0.00014(3) 


-0.00044(11) 


Ratio 


0.53(13) 


0.90(12) 


0.54(16) 


0.54(4) 
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FIG. 1. Difference with the exact results for the energy and the specific-heat in a 16^ lattice. We also plot Ai and A2. Pull 
circles correspond to Swendsen-Wang update, open ones to single-cluster and squares are from the Metropolis update. 
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FIG. 2. Simulation results for the energy, the specific- heat, Ai and A2 in a 16^ lattice. Dashed-lines for E and C„ are 

obtained from a minimization, excluding the SC-R250 data. Q is the probability of getting a larger value of x^- Full circles 
correspond to Swendsen-Wang update, open ones to single-cluster and squares are from the Metropolis update. 
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FIG. 3. Same as figure ^ for a 24'^ lattice. 
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